#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "omp.h" 
#include "mkl.h" 

#define IDF2(i, j, K) ((j)*(K) + (i))

double* getA16(int order)
{
	double *matrix = (double*)malloc(order*order*sizeof(double));
	for(int i = 0; i < order; ++i)
	{    
		for(int j = 0; j < order; ++j)
        {       
			matrix[IDF2(i,j,order)] = i>j? (order- i):(order - j);
        }
	}

	return matrix;
}

int main(int argc, char** argv)
{
	double* A = NULL;
	int i, j;
	int N = 6;
	if(argc >=  2)
	{
		N = atoi(argv[1]);
	}	
	A = getA16(N);

	int n = N;
	int info = -10;
	omp_set_num_threads(4); 

	printf("Start MKL fact \n");
	MKL_INT mkl_info;
	clock_t start_time, end_time;
	start_time = clock();
	dpotrf("L", &N, A,  &N, &mkl_info);
	end_time = clock();

	if(mkl_info)
	{
		printf("MKL FACTORIZATION FAILED mkl_info - %d\n", mkl_info);
		return EXIT_FAILURE;
	}

	//double *file_matrix = (double*)malloc(N*N*sizeof(double));
	//FILE * pFile;
 //   pFile = fopen ( "in/a16fact.dat" , "rb" );
	//fread (file_matrix , N*N , sizeof(double) , pFile );
	//fclose (pFile);


	//bool hasErrors = false;
	//for(int i = 0; i < N; i++)
	//{
	//	for(int j =0; j < N; j++)
	//	{
	//		if (fabs( file_matrix[IDF2(i, j, N)] - A[IDF2(i, j,N)] ) > 1e-6)
	//		{
	//			hasErrors = true;
	//			printf("For element [%i, %i] value of file %f value of mkl %f\n", i, j, file_matrix[IDF2(i, j, N)], A[IDF2(i, j,N)]); 
	//		}
	//	}
	//}

	//if(!hasErrors)
	//{
	//	printf("Comparing with MKL success\n");
	//}


	printf("\n Time : %.5f\n",(float)(end_time - start_time)/CLOCKS_PER_SEC);

	/*FILE * pFile;
	char file_name_buffer [50];
	sprintf(file_name_buffer, "in/a16_%d.dat", N);
    pFile = fopen ( file_name_buffer , "wb" );
	fwrite (A , N*N , sizeof(double) , pFile );
	fclose (pFile);*/

	free(A);
	//free(file_matrix);
}



